Import des packages

rm(list=ls())
library('dplyr')
library('tidyr')
library('VennDiagram')
library('grid')
library('futile.logger')
library("ggplot2")
library("gplots")
library('plotly')
library('tools')
library('cowplot')
library('UpSetR')

Import des VCF

Import des 3 vcf obtenus après lancement du pipeline (vcf_12x, vcf_24x, vcf40x) et import du vcf de référence (vcf_ref).

vcf_12x = read.table("/home/elodie/Documents/Analysis/gold_12x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_24x = read.table("/home/elodie/Documents/Analysis/gold_24x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_40x = read.table("/home/elodie/Documents/Analysis/gold_40x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_ref = read.table("/home/elodie/Documents/Elodie/Datas_ismael_ref/gold_on_data_elodie_ismael_decomposed_normalize_uniq.vcf")

Exploration des VCF et des variants

Nombre et type de variants par VCF

plot1 = ggplot(nombre_variants, aes(fill=type, y=vcf, x=nombre)) +
  geom_bar(position = "stack", stat = "identity") + 
  scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5)) + 
  labs(title = "Nombre et type de variant par VCF") 
ggplotly(plot1)

La distribution des types de variants semble équivalente pour chaque VCF. Cependant dans les 3 VCF générés par le pipeline on trouve des variants “Undetermined” mais pas dans le VCF de référence.

Que sont les Undetermined non retrouvés dans le VCF de référence ?

undetermined = vcf_all %>% filter(type == "undetermined")
print(undetermined)
##                                  variants chromosome position
## 1                       chr4_66388036_A_*       chr4 66388036
## 2 chr4_66440407_CATATATATACATATATATACAT_*       chr4 66440407
## 3                       chr4_66214417_A_*       chr4 66214417
## 4 chr4_66440407_CATATATATACATATATATACAT_*       chr4 66440407
## 5                       chr4_66473213_A_*       chr4 66473213
##                       ref alt                      format  Id  gene
## 1                       A   *      1/.:0,5:6:41:205,41,45 423 EPHA5
## 2 CATATATATACATATATATACAT   *  1/.:0,16:22:99:867,243,285 535 EPHA5
## 3                       A   *  1/.:0,11:15:99:717,152,110 110 EPHA5
## 4 CATATATATACATATATATACAT   * 1/.:0,22:31:99:1233,367,435 547 EPHA5
## 5                       A   *  1/.:0,4:15:99:1163,469,449 651 EPHA5
##           type     vcf
## 1 undetermined vcf_12x
## 2 undetermined vcf_24x
## 3 undetermined vcf_40x
## 4 undetermined vcf_40x
## 5 undetermined vcf_40x

Variants avec * en alt lié au vt décompose ?

pos66388036 = vcf_all %>% filter(variants == "chr4_66388034_GTA_G")
print(pos66388036)
##              variants chromosome position ref alt                        format
## 1 chr4_66388034_GTA_G       chr4 66388034 GTA   G         0/1:2,5:7:46:162,0,46
## 2 chr4_66388034_GTA_G       chr4 66388034 GTA   G       0/1:7,8:15:99:250,0,267
## 3 chr4_66388034_GTA_G       chr4 66388034 GTA   G     0/1:13,15:28:99:472,0,494
## 4 chr4_66388034_GTA_G       chr4 66388034 GTA   G 0/1:.:362:146,111:176,136:155
##    Id  gene     type     vcf
## 1 422 EPHA5 deletion vcf_12x
## 2 441 EPHA5 deletion vcf_24x
## 3 453 EPHA5 deletion vcf_40x
## 4 413 EPHA5 deletion vcf_ref
pos66388036bis = vcf_all %>% filter(variants == "chr4_66388036_A_G")
print(pos66388036bis)
##            variants chromosome position ref alt                   format  Id
## 1 chr4_66388036_A_G       chr4 66388036   A   G ./1:0,1:6:41:205,176,204 424
##    gene type     vcf
## 1 EPHA5  SNV vcf_12x

A cette position délétion de 2 bases (66388035+66388036) bien détectée dans les 4 VCF Substitution de A>G en 66388036 dans le VCF_12x (comme peu de reads, variant caller se trompe? Faux positif?)

pos66440407 = vcf_all %>% filter(position == 66440407)
print(pos66440407)
##                                  variants chromosome position
## 1 chr4_66440407_CATATATATACATATATATACAT_*       chr4 66440407
## 2 chr4_66440407_CATATATATACATATATATACAT_C       chr4 66440407
## 3 chr4_66440407_CATATATATACATATATATACAT_*       chr4 66440407
## 4 chr4_66440407_CATATATATACATATATATACAT_C       chr4 66440407
##                       ref alt                      format  Id  gene
## 1 CATATATATACATATATATACAT   *  1/.:0,16:22:99:867,243,285 535 EPHA5
## 2 CATATATATACATATATATACAT   C   ./1:0,6:22:99:867,631,654 536 EPHA5
## 3 CATATATATACATATATATACAT   * 1/.:0,22:31:99:1233,367,435 547 EPHA5
## 4 CATATATATACATATATATACAT   C  ./1:0,9:31:99:1233,875,889 548 EPHA5
##           type     vcf
## 1 undetermined vcf_24x
## 2     deletion vcf_24x
## 3 undetermined vcf_40x
## 4     deletion vcf_40x

Zone avec des reads ayant des dététions différentes qui se chevauchent dans les 4 BAM. Le VCF de référence n’appelle pas de variants dans cette zone. Zone difficile pour le variant caller, faire appel à un variant caller spécialiste des délétions?

pos66214417 = vcf_all %>% filter (position == 66214417)
print(pos66214417)
##            variants chromosome position ref alt                     format  Id
## 1 chr4_66214417_A_*       chr4 66214417   A   * 1/.:0,11:15:99:717,152,110 110
## 2 chr4_66214417_A_T       chr4 66214417   A   T  ./1:0,4:15:99:717,468,450 111
##    gene         type     vcf
## 1 EPHA5 undetermined vcf_40x
## 2 EPHA5          SNV vcf_40x
pos66214413 = vcf_all %>% filter (position == 66214413)
print(pos66214413)
##                variants chromosome position   ref alt                   format
## 1 chr4_66214413_TATAA_T       chr4 66214413 TATAA   T   0/1:1,9:10:15:375,0,15
## 2 chr4_66214413_TATAA_T       chr4 66214413 TATAA   T 0/1:4,11:15:99:450,0,110
##    Id  gene     type     vcf
## 1 108 EPHA5 deletion vcf_24x
## 2 109 EPHA5 deletion vcf_40x

Délétion de 4 bases de 66214413 à 66214417 vu dans les 4 BAM mais uniquement appelée dans VCF 24x et 40x. Non appelé dans le VCF de référence. Subtitution en 66214417 aussi un A>T vu dans les 4 BAM mais uniquement appelée dans le VCF 40x. Non appelé dans le VCF référence.

pos66473213 = vcf_all %>% filter (position == 66473213)
print(pos66473213)
##            variants chromosome position ref alt                      format  Id
## 1 chr4_66473213_A_C       chr4 66473213   A   C     1/1:0,18:18:54:720,54,0 635
## 2 chr4_66473213_A_*       chr4 66473213   A   *  1/.:0,4:15:99:1163,469,449 651
## 3 chr4_66473213_A_C       chr4 66473213   A   C ./1:0,11:15:99:1163,206,135 652
##    gene         type     vcf
## 1 EPHA5          SNV vcf_24x
## 2 EPHA5 undetermined vcf_40x
## 3 EPHA5          SNV vcf_40x

Dans le VCF de référence en position 66473209 à 66473213 : délétion de 4 bases appelée mais peu de reads avec la délétion et beaucoup plus de reads avec substitutions de A>C en 66473213. Substitutions appelée uniquement dans VCF 24x et 40x. Dans les zones où il y a des délétions et substitutions, le variant caller rencontre des difficultés?

Nombre de variants par gène et par VCF

plot2 = ggplot(nombre_variants_par_gene, aes(fill=gene, y=vcf, x=nombre)) +
  geom_bar(position="stack", stat="identity") +
  scale_fill_manual(values = c("#fdc086", "#91bfdb", "#b2df8a")) +
  labs(fill = "Nom des gènes", x = "Nombre de variants", y = "VCF") + 
  theme(panel.grid = element_blank(),plot.title = element_text(hjust = 0.5)) + 
  labs(title = "Nombre de variants par gène et par VCF") 
ggplotly(plot2)

Le nombre de variants est identique dans chaque VCF pour les gènes UTS2 et LPL mais le nombre varie pour le plus grand gène UPHA5.

Nombre et type de variants par gène et par VCF

nombre_variants_par_gene$vcf_nom_genes = paste(nombre_variants_par_gene$vcf, nombre_variants_par_gene$gene, sep="_")

plot3 = ggplot(nombre_variants_par_gene, aes(fill=type, y=vcf_nom_genes, x=nombre)) +
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
  labs(title = "Nombre de variants par gène, par VCF et par type de variant") +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") + 
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
  
ggplotly(plot3)

La répartition des types de variants par gène et par VCF est concordante. Aucune insertion dans UTS2 dans les 4 VCF.

Exploration des SNV

Nombre et type de SNV par VCF

plot4 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type), ], aes(fill = type, y = vcf, x = nombre))+
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  scale_fill_manual(values = c("#fc9272", "#9ebcda")) +
  labs(title = "Type de SNV par VCF", x = "Type de substitutions", y = "Nombre de substitutions") +
  labs(fill = "Type de SNV") +
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot4

En génétique humaine, il y a un environ 2 fois plus de transitions que de transversions. Sur le plot, on voit environ 2/3 de transitions et 1/3 de transversions ce qui est cohérent. Valeurs cohérentes dans tous les VCF.

Nombre et type de transitions/transversions par VCF et par gène

plot5 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type_snv), ], aes(fill = type_snv, y = vcf_gene, x = nombre)) +
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  labs(title = "Répartition des types de SNV par VCF et par gène", x = "Nombre", y = "VCF et gène") +
  scale_fill_manual(values = c("#a6cee3", "#1f78b4","#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a", "#ffff99", "#b15928")) +
  labs(fill = "Type de SNV", x = "Nombre", y = "VCF") +
  theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  theme(panel.grid = element_blank())
ggplotly(plot5)

Une répartition des types de SNV cohérente entre VCF pour chaque gène.

Exploration des délétions/insertions

plot6 = ggplot(vcf_all %>% filter (type == "insertion" | type == "deletion"), aes(x = taille_delins, y = vcf)) + 
  geom_boxplot(aes(fill=type)) + 
  theme_classic() + 
  labs(x = "Taille des delins", y = "VCF", title = "Boxplots de la taille des insertions et délétions par VCF") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Type de delins") +
  scale_fill_manual(values = c("deletion" = "#bebada", "insertion" ="#ffffb3"))
plot6

La distribution des tailles des délétions est équivalente entre les 4 VCF. Les VCF 12x, 24x et 40x contiennent plus de valeurs aberrantes pour les insertions et délétions que le VCF de référence. La taille maximale des insertions du VCF de référence est plus petite que pour les 3 VCF obtenus avec le pipeline. L’interquartile est plus petit également pour la taille des insertions du VCF de référence ce qui signifie que les valeurs des insertions sont plus regroupées autour de la médiane. Il y a moins de disparité des valeurs pour le VCF de référence.

plot7 = ggplot(vcf_all %>% filter (type == "insertion"), aes(x = taille_delins, y = vcf)) +
  geom_boxplot(aes(fill=gene)) +
  theme_classic() +
  labs(x = "Taille des insertions", y = "VCF", title = "Boxplots de la taille des insertions par VCF et par gène") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Noms des gènes") +
  scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))

plot8 = ggplot(vcf_all %>% filter (type == "deletion"), aes(x = taille_delins, y = vcf)) +
  geom_boxplot(aes(fill=gene)) +
  theme_classic() +
  labs(x = "Taille des délétions", y = "VCF", title = "Boxplots de la taille des délétions par VCF et par gène") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Noms des gènes") +
  scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))

plot_grid(plot7, plot8)

Exploration de la profondeur de lecture des variants

plot9 = ggplot(vcf_12x_24x_40x, aes(x = as.numeric(vcf_12x_24x_40x$DP), y = vcf)) + geom_boxplot(aes(fill=gene)) + theme_classic() + 
  theme_classic() +
  labs(x = "Profondeur de lecture des variants", y = "VCF", title = "Boxplots de la profondeur de lecture des variants par gène") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Noms des gènes") +
  scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086"))
plot9
## Warning: Use of `vcf_12x_24x_40x$DP` is discouraged.
## ℹ Use `DP` instead.

plot10 = ggplot(dp_12x, aes(x = DP, y = nombre)) +
  geom_histogram(stat = "identity", fill = "#a6cee3") +
  geom_vline(xintercept = mean(dp_12x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
  labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 12x") + 
  theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) + 
  theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#a6cee3"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot11 = ggplot(dp_24x, aes(x = DP, y = nombre)) +
  geom_histogram(stat = "identity", fill = "#1d91c0") +
  geom_vline(xintercept = mean(dp_24x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
  labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 24x") + 
  theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) + 
  theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#1d91c0"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot12 = ggplot(dp_40x, aes(x = DP, y = nombre)) +
  geom_histogram(stat = "identity", fill = "#253494") +
  geom_vline(xintercept = mean(dp_40x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
  labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 40x") + 
  theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5) ) + 
  theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#253494"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot_grid(plot10, plot11, plot12, ncol = 3, nrow = 1)

Exploration de la VAF des VCF 12x, 24x et 40x

On ne fera pas d’exploration de la VAF ni de la profondeur DP pour le VCF de référence car il n’y a pas d’explication aux valeurs de DP, AD, ADALL dans le VCF. Les valeurs ne sont pas cohérentes avec ce qui est observé dans le bam de référence associé au VCF dans IGV.

VAF minimale détectée dans chaque VCF

cat("VAF minimale dans VCF 12x :", min(vaf_12x$vaf), "%", "\n")
## VAF minimale dans VCF 12x : 16.7 %
cat("VAF minimale dans VCF 24x :", min(vaf_24x$vaf), "%", "\n")
## VAF minimale dans VCF 24x : 14.3 %
cat("VAF minimale dans VCF 40x :", min(vaf_40x$vaf), "%", "\n")
## VAF minimale dans VCF 40x : 16.3 %
plot13 = ggplot(vcf_vaf, aes(x = nombre_variants, y = cat_vaf , fill = vcf)) +
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  labs(title = "Répartition du nombre de variants par catégorie de VAF", x = "Nombre de variants", y = "Catégorie de VAF") +
  scale_fill_manual(values = c("#a6cee3", "#1d91c0", "#253494")) +
  labs(fill = "VCF") +
  theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  theme(panel.grid = element_blank()) +
  theme_classic()
ggplotly(plot13)

Exploration des variants en commun dans les 4 VCF

variants_list = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants), vcf_ref = unique(vcf_ref$variants))

plot14 = upset(fromList(variants_list), order.by = "freq", text.scale = 2,  point.size = 5, sets.bar.color = "#999999", main.bar.color = "#fc8d62")
plot14

plot15 = display_venn(variants_list, category.names = c("vcf_12x" , "vcf_24x" , "vcf_40x", "vcf_ref"), 
             fill = c("#a6cee3", "#1d91c0", "#253494", "#E69F00"), cat.cex = 1.2,
        cat.fontface = "bold",
        cat.default.pos = "outer",
        cat.dist = c(0.1, 0.1, 0.1, 0.1))

plot16 = ggplot(cat_vcf_all_nombre, aes(fill=type, x=nombre, y=categorie)) +
  geom_bar(position = "stack", stat="identity") +
  scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
   labs(title = "Nombre de variants par catégories et par type de variants") +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "Catégories de VCF") + 
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot17 = ggplot(cat_vcf_all_genes, aes(fill=gene, x=nombre, y=categorie)) +
  geom_bar(position = "stack", stat="identity") +
  scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086")) +
  labs(title = "Nombre de variants par catégories et par gènes") +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "Noms des gènes") + 
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))

plot_grid(plot16, plot17, ncol = 1, nrow = 2)

Exploration des variants en commun dans les 3 VCF générés par le pipeline

variants_list_12x_24x_40x = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants))

plot18 = display_venn(variants_list_12x_24x_40x, category.names = c("vcf_12x" , "vcf_40x", "vcf_24x"),
  fill = c("#a6cee3", "#1d91c0", "#253494"), cat.cex = 1.2,
        cat.fontface = "bold",
        cat.default.pos = "outer",
        cat.dist = c(0.1, 0.1, 0.1))

Exploration des métriques de qualité: datas générés par Hap.py d’Illumina

recall (sensibilité), précision (spécificité) et score F1. Ce sont des métriques qui permettent de voir si un modèle est performant. Le score F1 combien recall et précision Plus le recall est élevé, plus le modèle repère des positifs. Plus la précision est élevée, moins le modèle se trompe sur les positifs (moins il y a de faux positifs). Plus le score F1 est élevé, plus le modèle est performant.

Courbe ROC

happy_12x_summury = read.table("/home/elodie/Documents/Analysis/Happy_12x/Happy_12x/ref-12x.summary.csv", header=TRUE, sep=',')
happy_24x_summury = read.table("/home/elodie/Documents/Analysis/Happy_24x/Happy_24x/ref-24x.summary.csv", header=TRUE, sep=',')
happy_40x_summury = read.table("/home/elodie/Documents/Analysis/Happy_40x/Happy_40x/ref-40x.summary.csv", header=TRUE, sep=',')
plot_roc_snp = plot_data(roc_snp, is.PR=TRUE)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_roc_indel = plot_data(roc_indel, is.PR=TRUE)
plot_grid(plot_roc_snp, plot_roc_indel, nrow=2, ncol=1)